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Abstract 

This  research  investigated  the  effects  of  adding  a 
second-order  flux-difference-splitting  (FDS)  correction  term 
to  an  existing  computer  code  that  is  based  on  a  first-order 
FDS  algorithm.  It  was  determined  that  the  second-order 
algorithm  did  improve  the  accuracy  of  the  code  for  a  source 
flow  analysis,  but  second-order  behavior  could  not  be 
confirmed  by  the  error  convergence  patterns.  It  was  also 
discovered  that,  when  tested  across  an  oblique  shock  wave,  the 
second-order  correction  terms  had  minimal  influence  on  the 
accuracy  and  shock  capturing  ability  of  the  first-order 
accurate  FDS  method. 
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INVESTIGATION  OF  FLUX-DIFFERENCE-SPLITTING 


NUMERICAL  METHOD  IN  SUPERSONIC  NOZZLES 


I..  Introduction 


Purpose 

The  purpose  of  this  research  is  to  add  a  second-order 
accuracy  to  a  newly  developed  first-order  accurate  computer 
code  that  predicts  the  inviscid,  two-dimensional  (2D)  flow 
properties  within  and  around  the  nozzles  of  hypersonic 
aerospace  vehicles,  such  as  the  National  Aerospace 
Plane  (NASP).  The  code,  which  was  developed  by  Doty  (1),  is 
based  upon  an  innovative  differencing  technique,  known  as 
Flux-Difference-Splitting  (FDS).  The  flux-difference¬ 
splitting  numerical  method  was  created  by  Enquist  and 
Osher  (5:45-75).  Doty's  code  uses  a  first-order  accurate 
upwind  FDS  method,  operating  on  the  steady,  inviscid,  planar 
form  of  the  Euler  Equations  to  model  the  flow  field. 
Additionally,  the  fluid  is  modeled  with  a  perfect  gas 
assumption.  This  thesis  adds  second-order  correction  terms  to 
the  first-order  accurate  code  and  evaluates  the  results  by 
comparing  them  against  the  existing  first-order  accurate 
solution  and  known  exact  solutions.  The  exact  solutions  used 
for  this  comparison  are  supersonic  source  flow  and  shock  wave 
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reflection. 


Background 

Over  the  last  twenty  years,  computers  have  played  an 
increasing  role  in  the  design  of  aircraft.  Until  recently,  an 
aircraft  design  would  become  finalized,  and  then  a  scaled  down 
model  of  the  finished  design  would  be  built  and  tested  in  a 
wind  tunnel.  This  approach  proved  to  be  very  expensive,  time 
consuming,  and  inefficient.  Additionally,  it  lacked  the 
flexibility  needed  to  improve  the  design  as  contract  or 
performance  requirements  changed,  or  to  allow  optimization  of 
the  design  during  the  design  phase.  Now,  with  the  advancement 
of  supercomputers,  engineers  have  the  capability  to 
mathematically  model  ai.d  simulate  flights  of  the  current 
unfinished  design,  making  it  possible  to  optimize  the  design 
before  it  is  finalized. 

According  to  Barthelemy  (2:6-9),  the  NASP  will  require 
technology  that  is  a  quantum  leap  ahead  of  the  technology  that 
is  used  in  today's  aircraft  and  spacecraft,  and  will  expend 
enormous  amounts  of  government  and  contractor  rt-ourcas.  A 
sketch  of  a  typical  NASP  type  vehicle  can  be  seen  in  Figure  1. 
This  type  of  vehicle  will  be  powered  by  a  supersonic 
combustion  ramjet  (SCRAMJET)  and,  theoretically,  will  be  able 
to  cruise  at  speeds  of  up  to  Mach  25. 

The  computer  code  was  written  as  a  "user  friendly"  tool 
that  design  engineers  could  utilize  to  predict  flow  parameters 
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in  a  NASP  type  nozzle  to  enhance  design  optimization.  These 
predicted  flow  parameters  are  shock  wave  and  slip  stream 
locations,  pressures,  temperatures,  and  velocities  throughout 
the  plane  of  the  nozzle. 

Figure  2  shows  the  general  configuration  of  a  NASP  type 
nozzle.  In  Figure  2,  region  1  is  the  combustor  exit. 
Region  2  is  the  external  air  flow  that  passes  under  the 
engines.  The  exhaust  flow,  in  region  3a,  and  the  external  air 
flow,  in  region  3b,  are  separated  by  a  contact  surface  that 
originates  from  the  engine  cowl,  as  illustrated  in  the  figure. 
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Typical  hypersonic  vehicle  (3:4). 
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Region  (i) 
External  flow 


Contact  surface  -7 


Region  © 

External  flow  interaction 


y 


Fictitious  lower  boundary 


X 

Figure  2.  Expanded  view  of  nozzle  and  cowl  section  (3:5). 
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II .  Governing  Equations  and  Coordinate  Transformation 


Governing  Equations 

According  to  Anderson  et  al.  (1:235-236),  for  flows  with 
sufficiently  high  Reynolds  numbers  the  viscous  and  heat 
transfer  effects  are  confined  to  a  thin  boundary  layer  near 
the  wall.  For  a  large  ducted  supersonic  flow,  such  as  flow  in 
a  NASP  type  nozzle,  these  effects  can  be  neglected  for 
preliminary  analysis  and  design.  With  these  two  assumptions 
applied  to  the  full  Navier-Stokes  equations,  and  neglecting 
body  forces,  the  governing  equations  for  this  2D,  steady, 
compressible,  inviscid,  adiabatic  fluid  flow  can  be  reduced  to 
the  inviscid  Euler  Equations: 


dx  dy 


(1) 


where 


S  = 


pu 

pu^  +  P 
puv 

u(pe  +  P) 


pv 

pvu 

pv^  +  P 
v(pe  +  P) 


(2) 


These  equations  are  derived  from  the  conservation  laws: 
conservation  of  mass,  conservation  of  momentum,  and 
conservation  of  energy. 
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Coordinate  Transformation 


Equation  (1)  applies  to  the  physical  domain  and  can  be 
transformed  into  the  computational  domain,  where  it  is  more 
convenient  to  numerically  solve.  Figure  3  shows  both  the 
physical  and  computational  coordinate  systems  for  the  nozzle. 
To  accomplish  this  transformation,  the  following  mapping  is 
applied: 


C  =  =  T)  (x,y) 


(3) 


Here,  the  axial  coordinate,  x,  maps  directly  into  the 
computational  domain,  but  the  normal  coordinate,  y,  is 
transformed  by  a  nonlinear  function.  To  determine  what  these 
functions  are,  the  chain  rule  of  multivariable  calculus  must 
be  applied: 


ao  _  dc  ao  ^  ao  _  ,  ao  ^  ao 
ax  ax  ac  ax  ac  ari 

ILL  =  .  an  ao  _  .  ao  ^  ao 

dy  dy  ac  ay  a^  ac  an 


From  Eg  (3),  it  can  be  seen  that: 
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c.  =  1 


Cy  =  0 


(6) 


Combining  Eqs  (4),  (5),  and  (6)  gives  the  chain  rule  form  that 
is  used  to  transform  the  governing  equations  from  the  physical 
to  the  computational  domain: 


ao  _  ao  ^  „  ao 
dx  ac  an 


(7) 


ao  _  80 

dy  dn 


(8) 


Applying  Eqs  (7)  and  (8)  to  the  original  governing 
equations,  Eqs  (1)  and  (2),  gives  the  transformed  governing 
equations  in  the  computational  domain: 


d{B)  ^  ^  d{E) 


+  V 

Sy 


d{F) 

an 


=  0 


(9) 


Equation  (9)  can  be  written  in  a  more  convenient  form  as 


d(B) _ d(S)  diF) 

ac  an  an 


(10) 


For  further  information  on  the  governing  equations  or  the 
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transformation  of  the  equations,  see  Appendix  G  in 
reference  3. 
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Ill . 


FLUX-DIFFERENCE-SPLITTING 


Riemann  Problem 

The  Riemann  problem  is  the  heart  of  the  FDS  method.  The 
flowfield  must  be  discretized  before  the  Riemann  problem  is 
applied.  The  first  step  in  the  process  is  to  model  the 
general  flow  property,  i|f  in  Figure  4a,  which  has  an  arbitrary 
spatial  distribution,  as  a  series  of  nodes.  Each  node 
represents  a  localized  region  of  uniform  flow.  The  Riemann 
problem  assumes  that  these  localized  uniform  flow  regions,  at 
j  and  j+1  in  Figure  4a,  extend  to  a  distance  half-way  between 
the  two  nodes,  and  a  discontinuity  is  assumed  to  exist  at  the 
midpoint,  j+1/2. 

Waves  are  generated  by  the  discontinuity  at  j  +  1/2.  These 
waves  can  be  seen  in  Figure  4b.  Waves  (1)  and  (3)  can  either 
be  shock  waves,  expansion  waves,  or  a  combination  of  one  of 
each  type  of  wave  depending  on  the  local  flow  configuration. 
Wave  (2)  is  a  contact  surface.  These  three  waves  separate  the 
flow  into  four  regions  of  uniform  flow,  with  regions  6  and  0 
having  values  equal  to  known  values  at  nodes  above  and  below 
the  Riemann  location,  respectively.  Regions  2  and  4 
correspond  to  unknown  values  at  the  Riemann  location  that  must 
be  solved  for. 

Three  different  methods  exist  to  solve  the  Riemann 
problem  (3:10-13).  The  first  method  is  an  exact  solution  to 
the  Riemann  problem  where  an  iteration  is  performed  using  the 
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oblique  shock  wave  relationships  across  shock  waves,  and 
Prandtl-Meyer  expansion  relationships  across  expansion  waves. 
The  second  method  approximates  the  shock  wave  as  an  isentropic 
compression  and  again  requires  iteration  using  Prandtl-Meyer 
relationships  across  both  expansions  and  compressions.  The 
third  method  is  a  linearized  approximate  solution  that  treats 
the  local  Riemann  problem  as  isentropic.  This  method  uses  a 
linearized  form  of  the  Prandtl-Meyer  relationships  to 
analytically  solve  across  both  expansions  and  compressions. 
Only  the  exact  solution  is  used  in  this  research. 

Exact  Solution  to  the  Riemann  Problem.  The  Riemann 
problem  is  solved  exactly  for  the  properties  in  regions  2  and 
4,  as  shown  in  Figure  4b,  by  using  oblique  shock  and  Prandtl- 
Meyer  expansion  relationships.  This  is  an  iterative  process 
in  which  the  parameters  in  regions  2  and  4  are  calculated 
separately,  but  there  are  two  known  conditions.  First,  the 
pressure  in  region  4  must  match  that  of  region  2  since  the 
contact  surface  cannot  support  a  normal  pressure  gradient. 
The  second  known  condition  is  that  the  flow  angle  in  region  4 
is  equal  to  the  flow  angle  in  region  2  with  the  two  angles 
equal  to  that  of  the  contact  surface. 

Shock  Wave  Relationships.  For  the  example  in 
Figure  4,  the  pressure  at  j+1  is  greater  than  the  pressure  at 
node  j.  Therefore,  wave  1  is  a  shock  wave,  and  must  be  solved 
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using  the  oblique  shock  wave  relationships.  The  following 
shock  wave  relationships  are  written  with  0  corresponding  to 
known  values  upstream  in  Riemann  region  0,  and  2  corresponding 
to  values  that  must  be  solved  for  downstream  in  Riemann 
region  2.  In  the  oblique  shock  example  shown  in  Figures  4  and 
5,  the  pressure  in  region  2  is  greater  than  the  pressure  in 
region  0.  The  following  oblique  shock  wave  equations  are 
taken  from  Zucrow  (9:359-360). 

5^  =  0a  -  00  (11) 

1 

tan(62) 

P  =  c  -  62  (13) 


Y+l 


Wn  sin^(e)  -  1 


-  1 


tan(e) 


(12) 


—  =  — Mq  sin^(e) 
Po  Y  +  1 


Y  -  1 

Y  +  1 


(14) 


£2  ^  tan(e)  ^  (y  ^  1)M^  sin^(e) 

Po  tan(P)  2  +  (y  -  1)Mq  sin^(c) 
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2 


(16) 


Yi 

^0 


sin (e) 
sin(P) 


+  DMo  sinMe) 


.  Y__i 

Y  +  1 


tan (e)  _  2  / _ i  ^  y  ~  i’) 

tan(p)  Y  ^  1  [m^  sin2(P)  2  j 


Uj  =  cos  (62) 


(18) 


V2  =  V2  sin (62) 


(19) 


Expansion  Wave  Relationships.  For  the  example,  we 
are  assuming  isentropic  expansion,  and  using  Prandtl-Meyer 
expansion  relationships.  In  this  analysis,  the  subscripts 
correspond  to  the  applicable  Riemann  region.  See  Figures  6 
and  7  for  a  graphical  representation  of  the  example.  The 
following  equations  are  used  across  an  expansion  fan: 

64  =  04  -  06  (20) 


b  = 


Y  ^  1 

Y  -  1 


1/2 


(21) 


14 


=  b  arc  tan  -  arctan^A^^  -l  (22) 


V4  =  Vg 


-  64  =  arctan^y/M^  -1  -  arctan^W^  -1  (23) 


■  4  _ 


.  V--  .1 


1  + 


(Y  -  1) 


(24) 


fi  = 
Pe 


(Y  -  1) 


(25) 


= 


Y  P4 


1/2 


(26) 


U4  =  W4  003(64) 


(27) 


V4  =  34  sin  (64) 


(28) 


Remember,  as  stated  above,  in  Riemann  regions  2  and  4  the  flow 
angles  and  pressures  must  match: 
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(29) 


04  =  02 


^4 


=  Po 


(30) 


This  is  an  iterative  process,  and  the  flow  turning  angles 
may  change  from  compression  to  expansion,  expansion  to 
compression,  or  experience  no  turning.  Therefore,  the  flow 
turning  angle  should  be  checked  after  each  iteration  to 
determine  the  analytical  mode  that  will  be  used  next.  The 
iteration  process  is  as  follows  (3:166-167): 

1.  The  initial  conditions  for  regions  6  and  0  are  known  and 
remain  fixed. 

2.  The  flow  angle  for  wave  (2)  is  guessed.  (The  average  of 
the  flow  angles  between  regions  6  and  0  is  usually 
sufficient  as  a  first  guess). 

3.  Solve  the  shock  wave  problem  for  the  pressure  in  region  2. 

a.  The  flow  deflection  angle  is  calculated  as  the 
change  in  the  flow  angle  between  regions  0  and  2  for  the 
shock  from  Eq  (11). 

b.  For  the  assumed  turning  angle,  the  oblique  shock  wave 
angle  is  iteratively  calculated  from  Eq  (12). 

c.  The  static  pressure  in  region  2  is  then  calculated 
from  Eq  ( 14 ) . 
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4.  Solve  the  expansion  wave  problem  from  regions  6  to  4 . 

a.  Using  Eg  (29)  the  flow  angle  in  region  4  is 
reguired  to  be  the  same  as  the  flow  angle  in  region  2. 

b.  The  flow  deflection  angle  is  then  calculated  as  the 
change  in  the  flow  angle  between  regions  6  and  4  for  the 
expansion  wave  using  Eg  (20). 

c.  The  Prandtl-Meyer  angle  is  calculated  for  region  6 
using  Eg  (22)  and  remains  fixed. 

d.  The  Prandtl-Meyer  angle  is  calculated  for  region  4 
using  the  first  of  Eqs  (23). 

e.  The  Mach  number  in  region  4  is  iteratively  calculated 
from  the  second  of  Eqs  (23). 

f.  The  static  pressure  in  region  4  is  calculated  from 
Eg  (24). 

5.  Check  for  consistent  solutions. 

a.  The  static  pressure  in  regions  4  and  2  must  match 
across  the  contact  surface  according  to  Eg  (30). 

b.  If  the  pressures  are  equal  (or  within  a  specified 
tolerance),  the  iteration  is  complete. 

c.  If  the  pressures  are  not  equal,  a  new  flow  angle  is 
guessed  and  the  iteration  continues  from  step  3. 

6.  Calculate  remaining  properties  from  the  shock  and 
expansion  wave  relations. 
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Flux  Differencing 


Before  the  fluxes  can  be  differenced,  they  must  be 
computed.  This  is  accomplished  by  using  the  primitive 
variables  that  were  calculated  by  the  exact  solution  to  the 
Riemann  problem.  Recall  that  these  primitive  variables  were 
computed  at  the  Riemann  locations,  which  are  positioned  half¬ 
way  between  the  nodes,  as  shown  in  Figure  8.  The  Riemann 
fluxes  must  be  calculated  for  all  four  components  of  the  E  and 
F  vectors,  in  each  of  the  four  Riemann  regions,  and  are 
recombined  as  defined  by  Eqs  (1)  and  (2),  repeated  here  for 
convenience ; 


i?  .  ^  =  0 

dx  by 


(1) 


where 


pu 

pu^  +  P 
puv 
[u(pe  +  P) 


F  - 


pv 

pvxi 

pv^  +  p 

v(pe  +  P) 


(2) 


As  an  example,  the  first  E  flux  component,  El,  would  be 
computed  at  each  Riemann  location  for  all  four  Riemann 
regions : 


18 


(El)o  =  PoUq 


(31) 


(ED^^p^u^  (32) 

(El),  =  p,U4  (33) 

(Ei),  =  PgUg  (34) 

After  the  fluxes  have  been  calculated  at  all  the  Riemann 
locations  on  the  vertical  plane,  see  Figure  8,  they  must  be 
locally  differenced  across  waves  (1),  (2),  and  (3). 

Differencing  El  across  wave  (1)  corresponds  to  differencing 
the  density/axial  velocity  product  between  regions  0  and  2. 
This  can  be  shown  mathematically  as: 

=  (El).  -  (EI)o  =  92^2  -  Po^o  (35) 

Likewise,  differencing  El  across  waves  (2)  and  (3): 

idEl)„^^^2  =  (El),  -  (El)^  =  p,U,  -  P2U2  (36) 
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{dBl) 


wave  3 


( El )  g  (  El )  4  P  6  ^6  P  4  ^4 


(37) 


Summing  the  flux  differences  across  all  three  waves  gives  the 
total  flux  difference  at  the  Riemann  location,  j+1/2: 

(dEl)  -  [  (dEl)  ^  +  (dEl)„^y^2  wave  I  ^  jf*l/2 


The  differences  of  the  other  three  E  components  and  the  four 
F  components  are  computed  in  a  similar  fashion. 

At  this  point,  the  FDS  approach  is  still  similar  to  the 
finite-difference  approach  in  that  Eg  (38)  represents  the 
total  difference  between  nodes  j  and  j+1.  This  is  where  the 
FDS  method  differs  from  the  finite-difference  method  in  that 
for  FDS  the  differences  are  split  into  positively  and 
negatively  biased  components. 

Splitting  the  Flux  Differences 

Splitting  the  flux  differences  is  a  directional  biasing 
of  the  flux  differences.  As  a  result,  only  relevant 
information  is  received  at  each  downstream  node,  as 
illustrated  in  Figure  9. 

Splitting  the  flux  differences  is  accomplished  by 
breaking  up  the  flux  differences,  that  were  calculated  at  each 
Riemann  location,  into  positive  and  negative  components.  For 
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the  El  flux  difference,  splitting  across  each  wave  would  give: 


(dEl),,,,.  =  +  {dEl-)„,,,^  (39) 

(dEl),,,,.  =  ^  (dEl-)„,,,,  (40) 

(dEl)„,,,,  =  {dEl*),,,,,  ^  {dEl-)^,,,  (41) 


Here,  the  "+"  superscript  indicates  that  this  quantity  is  the 
portion  of  the  flux  difference  that  passes  information  in  the 
positive  y  direction.  Inversely,  the  superscript  only 

passes  information  in  the  negative  y  direction.  As  before, 
the  split  flux  differences  can  be  added  over  all  three  waves 
to  get  the  total  positive  and  negative  flux  differences  at  the 
Riemann  location  for  all  E  and  F  components. 

(  dEl )  =  ((dEl)  1  +  ( d£'i  )  lave  2  +  (dEl)  ^  ^ 

idEl)  jt.i/2  ~  {idEl)  1  (dEl)  2  *  (dEl)  (43) 


Again,  this  procedure  would  be  repeated  for  the  other  three  E 
and  four  F  vector  components. 
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Values  are  known  at  computational  plane  i,  and  the  split 
flux  differences  are  used  to  pass  information  to  the  next 
plane,  i+1.  See  Figure  9  for  a  graphical  representation.  The 
information  is  passed  in  the  direction  of  characteristics  and 
streamlines  in  the  computational  domain.  These  characteristic 
slopes  are  given  by  (dy/dx)  in  the  physical  domain,  but  in  the 
computational  domain  they  are: 


^  dC  dx  dx  dx  dy  dx 


(44) 


^2  “  'lx 


dx 


Aj  =  n, 


dy] 

dx 


(45) 


(46) 


where 


dyl  _  uv  -  a^yjM^-1 
dxJi  u2  _  ^2 


(47) 
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l2 


V 

U 


(48) 


'  dy 
dx 


'  dy}  _  uv  +  -  1 

.  dx\^  u2  -  a2 


(49) 


Equations  (47),  (48),  and  (49)  can  be  manipulated  to  give  the 
characteristic  slopes  in  the  physical  domain  in  terms  of  flow 
angle  and  Mach  angle,  see  Doty  (4): 


=  tan(0  -  a) 


(50) 


^1  =  tan  (6) 
dx\2 


(51) 


=  tan(0  +  a) 


(52) 


Maximum  Step  Size  for  Marching  Algorithm 

The  step  size,  for  this  steady  problem,  corresponds  to 
the  largest  axial  distance,  from  plane  i  to  plane  i+1,  that 
can  be  made  while  keeping  the  solution  stable,  and  must  be 
determined  at  each  computational  plane.  According  to  Anderson 
et  al.  (1:76),  the  CFL  condition  requires  that  the  analytic 
domain  of  influence  must  lie  within  the  numerical  domain  of 
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influence,  as  illustrated  in  Figure  10.  Therefore,  any  step 
size  larger  than  this  amount  would  cause  instability.  After 
all  the  characteristic  slopes  and  stream  line  locations  have 
been  computed  at  the  computational  plane,  i,  a  local  step 
size,  AC/  can  be  determined  for  each  nodal  location,  j.  The 
local  step  size  is  egual  to  the  minimum  axial  distance  of  the 
intersection  of  the  streamline  from  node  j  and  either  the 
negatively  biased  characteristic  at  j+1  or  the  positively 
biased  characteristic  from  j-1.  The  magnitude  of  AC  for  the 
computational  plane,  i,  is  equal  to  the  minimum  local  AC  on 
that  plane.  Also,  for  this  research,  a  Courant  number 
multiplier  of  0.99  was  used. 
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Figure  4,  Riemann  problem  for  planar  supersonic  flow  (3:14). 
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Figure  5.  Oblique  shock  wave  geometry  (3:178). 
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(a).  Pressure  distribution.  (b)-  Wave  pattern. 

Figure  7.  Specific  Riemann  problem  (3:179) 


Figure  10.  Step  size  determination  (3:157). 
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IV.  First-Order  Accurate  Upwind  Flux-Dif f erence-Split 

Numerical  Algorithm 

First-Order  Accurate  Interior  Point  FDS  Approximations 

Recall  that  the  governing  equations  for  steady,  inviscid, 
adiabatic,  planar  flow  in  computational  space  are  given  by: 


d(S)  _  ^  d(S)  _  ^  d{F) 

dc  ~  dr\  dr\ 


(10) 


These  equations  are  operated  on  by  a  flux-difference-splitting 
in  the  normal  direction,  and  finite  differencing  in  the  axial 
direction.  The  governing  equations  can  be  rewritten  into 
their  approximate  computational  form,  shown  in  Eq  (53), 
(3:190). 


^AE)  ^AE)  ^AE) 


(S3) 


Here,  for  computational  convenience. 

At,  =  1  (54) 

Substituting  Eq  (54)  into  Eq  (53),  and  rearranging  gives: 
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Ai(^)  =  -AC  Tlx  A^(S)  -AC  Tiy  iF) 


(55) 


The  "i"  subscript  indicates  a  finite-difference  operator  in 
the  axial  direction.  In  this  case,  a  two-point,  upwind, 
first-order  accurate,  finite-difference  operator  is  applied 
between  planes  i  and  i+1.  The  "j"  subscript,  in  Eq  (55), 
signifies  that  a  FDS  operator  is  applied  in  the  normal 
direction.  Here,  a  two-point,  upwind,  first-order  accurate, 
FDS  operator  is  applied  in  both  the  positive  and  negative  y 
directions  at  node  j  on  plane  i,  resulting  in  a  three  point 
stencil  in  the  normal  direction,  as  shown  in  Figure  11,  In 
Eq  (55),  the  finite-difference  operator  is; 

A^IE)  =  -  Ei  (56) 


The  flux-difference  operators  are: 


AjiE)  =  [c®;.i/a  +  dE/.i/a] 


(57) 


and 


A,.(F)  =  [dF/.i/a  +  dF/.i/a] 


(58) 


Note  how  only  the  positive  biased  information  from  the  Riemann 
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location  below  j,  at  j-1/2,  and  the  negative  biased 
information  from  the  Riemann  location  above,  at  j+1/2,  is  used 
to  determine  the  new  values  at  the  next  computational  plane, 
i+1.  Substituting  Eqs  (56),  (57),  and  (58)  into  Eg  (55)  and 
rearranging  gives: 


-  AC  Tiy  [cirT-i/a  ^ 


(59) 


Equation  (59)  represents  the  first-order  differencing  equation 
in  chain  rule  conservation  form.  The  chain  rule  conservation 
form  uses  metrics  that  are  calculated  at  the  nodes.  The  weak 
conservation  law  form,  which  is  another  method  for  solving  the 
governing  equations,  can  be  found  in  Appendix  A. 


First-Order  Accurate  Boundary  Point  FDS  Approximations 

When  calculating  values  at  the  upper  boundary,  of 
plane  i+1,  there  is  no  Riemann  location  above  the  boundary  to 
use,  as  can  be  seen  in  Figure  12.  Therefore,  only  physical 
information  at  Riemann  location  j-1/2  can  be  used.  Thus,  the 
split  flux  differences  coming  from  j+1/2  do  not  exist  and  are 
deleted  from  Eg  (59)  for  the  upper  boundary.  For  the  upper 
boundary.  Eg  (59)  becomes: 
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=  */  -  -  AC  ny(d^*>j-l/2 


(60) 


For  the  lower  boundary: 

if*  =  1/  -  A{  t|,  (dilj.irt  -  AC  1),  <«A> 

After  the  partial  solution  at  plane  i+1  has  been  computed 
using  Eqs  (60)  and  (61),  the  conservative  variables  are 
decoded  into  the  primitive  variables  p,  u,  v,  p,  and  pe.  The 
resultant  decoded  solution  does  not  necessarily  obey  the 
inviscid  velocity  tangency  condition  at  the  wall.  Therefore, 
to  be  physically  consistent  with  the  geometry,  the  primitive 
variables  need  to  be  corrected,  and  a  wave  corrector  is 
applied  to  the  solution  to  turn  the  flow  parallel  to  the  wall. 
To  accomplish  this  turning,  a  new  deflection  angle,  5,  is 
specified  as: 


®  ®soi  ®i,aii 

where 

is  the  computed  flow  angle 
^waii  wall  angle 


(62) 


After  the  deflection  angle  is  known,  the  flow  can  be  turned 
parallel  to  the  wall  by  one  of  two  methods.  If  a  compression 
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wave  is  required  to  turn  the  flow,  the  oblique  shock  wave 
relationships  are  used.  A  shock  wave  angle,  e,  can  be 
iterated  using  Eq  (12),  and  the  corrected  pressure,  density, 
and  velocity  components  are  computed.  If  and  expansion  is 
required  to  turn  the  flow,  the  Prandtl-Meyer  expansion  wave 
relationships  are  used.  Equation  (23)  would  be  iterated  to 
determine  the  downstream  Mach  number.  Once  the  Mach  number  is 
known,  the  primitive  variables  can  be  computed. 
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i+1 

X  or  ^ 


Figure  11.  Stencil  for  first-order  accurate  FDS  upwind 
method  (3:196). 


Figure  12.  Stencil  for  first-order  accurate  upper  solid 
wall  boundary  point  (3:196). 
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V.  Second-Order  Accurate  Center-Spaced  Flux-Difference- 

Split  Numerical  Algorithm 


Introduction 

The  first-order  accurate  method  is  the  basis  upon  which 
the  second-order  accurate  center-spaced  method  is  built. 
Recall  that  the  solution  is  advanced  with  first-order  accuracy 
using  Eq  ( 59 ) : 


if"  =  i/  -  AC 

-  AC  Hy  [d^j  -i/a  ciFj+i/a] 


(59) 


Equation  (59)  can  be  rewritten  as 

(63) 

where 

^  -  AC  Tlx  ((d»%-l/2  (C»')>l/2] 

(i)  ^  -  AC  Tiy  [(d''%-i/2  idr-)j.u2] 


In  the  above  equations,  the  Ixy  subscript  indicates  that  these 
quantities  are  first-order  accurate  corrections  in  the  x  and 
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Y  directions,  and  are  functions  of  the  sums  of  the  flux 
differences . 


Second-Order  Accurate  Center- Spaced  Interior  Point 
FDS  Approximations 

Second-order  accuracy  can  be  implemented  by  two  different 
methods.  First,  and  most  obvious,  is  to  use  a  second-order 
accurate  differencing  method  for  the  derivative  approximations 
in  the  governing  equations.  The  second  method,  which  is  used 
here,  uses  a  first-order  accurate  solution  and  adds  a  second- 
order  corrector  to  it.  Thus,  the  second-order  accurate 
solution  is  comprised  of  the  first-order  accurate  solution  and 
second-order  corrections. 

Doty  has  shown  that  a  second-order  correction  applied  to 
the  linear  hyperbolic  convection  equation,  which  is  often  used 
to  model  the  essential  characteristics  of  the  nonlinear  system 
of  equations,  has  a  modified  equation  that  demonstrates 
second-order  accuracy  (3:108).  For  a  complete  description  of 
this  process,  refer  to  Appendices  A,  B,  and  L  in  reference  3. 

The  second-order  accurate  equation  can  be  written  in 
shorthand  notation  as; 

=  s/  +  ^  (S) ,  +  A  •  (Jff) ,  +  A  .  (Jff) , 

"j  ixy  “7 '**' 2y  “7 2x  (66) 

+  Aj(^2y  ^7^^2x 
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The  second-order  normal  corrections  are: 


^  -  I  A{  1).  [(d«*>i.i/2  -  (caij-i/i] 

*  I  AC  t|, 

and 

■  -  I  AC  n,  [(d»’*>j.i/2  - 

*  I  AC  n,  [(c!»-)j.„a  -  (df-)j.u,] 


The  second-order  axial  corrections  are: 


Aj  <*)  =  ■§  AC^  11,  [((A  -  (A  ■!«*), -1/2). 

1-  ((A  d**)...,.  -  (A 
t  ((A  -  (A 

-  j  AC“  t|,  [((|A|  d*-).,./.  -  (|A|  d*-).-..,.), 

»  ((|A|  d*-).,./.  -  (|A|  d«-).../.)j 
*  {(|A|  d*-).,.„  -  (|A|  dB-)_,../.),] 


=  -i  A{^  n,  [((A  -  (A  dFVi,.). 

+  ((A  dP*),-.,/,  -  (A  dP*)j-i/i); 

♦  ((A  dP-)^.„j  -  (A  dPT-irt),] 

^  (70) 

-  I  A{^  i,^  [((|A|  dP-)j„„  -  (|A|  dr-)j.,/,)^ 

•  ((|A|  dr-)j.,„  -  (|A|  dp-)j_i,j)^ 
t  ((|A|  dP-)j,i„  -  (|A|  dPIj-irt),] 


Note  that  the  second-order  corrections  are  differences  of  the 
split  flux  differences,  and  both  the  axial  and  normal 
correction  are  required  to  achieve  overall  second-order 
accuracy.  As  can  be  seen  in  Figure  13,  the  second-order 
corrections  do  not  pass  information  in  the  same  manner  as  the 
first-order  correction.  The  second-order  method  passes 
positive  information  from  node  j+1/2,  and  negative  information 
from  node  j-1/2.  This  violates  the  rules  that  only  positive 
biased  information  is  used  from  the  Riemann  location  below 
node  j ,  and  only  negative  biased  information  is  used  from  the 
Riemann  location  above  node  j •  This  is  not  a  problem  due  to 
the  small  magnitude,  and  influence,  of  second-order  correction 
terms,  Pandolfi  (7:606). 
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Second-Order  Accurate  Center-Spaced  Boundary  Point 
FDS  Approximations 


The  second-order  boundary  correction  terms  are  computed 
in  a  manner  similar  to  the  first-order  boundary  solution.  The 
difference  is  that  information  is  required  at  a  fictitious 
Riemann  location  just  outside  the  boundary,  at  j+1/2,  for  the 
upper  boundary  in  Figure  13.  This  fictitious  location  is 
assigned  information  so  that  a  difference  of  a  split  flux 
difference  can  be  computed  for  the  Riemann  location  just 
inside  the  boundary,  at  j-1/2.  As  an  example  on  the  upper 
boundary,  to  get  a  difference  of  a  split  flux  difference  at 
j-1/2,  the  split  flux  differences  at  locations  j-1/2  and  j-3/2 
must  be  extrapolated  to  j+1/2.  According  to  Pandolfi  (7:607), 


a  linear  extrapolation  is  sufficient: 

.  2  (ds-)j.i/,  -  (71) 

=  2  (72) 

On  the  lower  boundary: 

J-1/2  '  2  (d'-lj.i/j  -  (<«'■)  j.i/2  (’*> 
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As  an  example  on  the  upper  boundary,  after  the  split  flux 
differences  are  determined  for  the  fictitious  Riemann 
location,  j+1/2,  a  difference  of  a  split  flux  difference  can 
be  computed  at  j-1/2.  The  second-order  correction  is  computed 
in  a  procedure  similar  to  that  of  the  first-order  solution, 
except  instead  of  using  split  flux  differences  at  j-1/2,  the 
second-order  method  uses  differences  of  split  flux  differences 
at  j-1/2  to  compute  the  second-order  corrections.  Once  the 
second-order  correction  is  known,  it  is  added  to  the  first- 
order  solution.  The  same  method  is  repeated  on  the  lower 
boundary.  See  Chapter  IV. 
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Figure  13.  Stencil  for  second-order  accurate 
center-spaced  FDS  method  (3:213). 
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Figure  14.  Stencil  for  second-order  accurate 
center-spaced  FDS  method  at  upper  boundary  (3:214). 


41 


VI .  Results  and  Discussion 


Two  different  investigations  were  performed  to  determine 
the  effects  of  the  second-order  correction  terms.  The  first 
investigation  was  a  grid  refinement  study  using  a  supersonic 
source  flow.  The  second  was  a  method  comparison  using  an 
oblique  shock  wave  study. 

Grid  Refinement  Study 

Before  the  grid  refinement  study  is  discussed,  a 
convention  must  be  defined.  In  this  discussion,  methodXX  will 
be  used  to  describe  a  differencing  scheme  used  in  the 
analysis,  where  the  first  X  corresponds  to  the  differencing 
scheme  used  at  an  interior  location,  and  the  second  X  would 
indicate  the  method  used  on  the  boundaries.  The  definition  is 
as  follows; 

1-  corresponds  to  a  first-order  upwind  differencing 

2-  corresponds  to  a  second-order  central  differencing 
Or,  stated  directly: 

raethodll-  1st  order  interior,  1st  order  boundary 
method21-  2nd  order  interior,  1st  order  boundary 
raethod22-  2nd  order  interior,  2nd  order  boundary 
These  conventions  will  be  used  from  this  point  forward. 

The  grid  refinement  study  was  accomplished  with  a  source 
flow  comparison,  using  the  initial  conditions  in  Table  I. 
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Table  I.  Initial  conditions  for  source  flow  analysis. 


Property 

Mach=1.01 

static  pressure  (Pa) 

101,325 

static  temperature  (K) 

298.0 

specific  heat  ratio 

1.4 

gas  constant  (J/Kg/K) 

287.0 

An  example  of  the  source  flow  geometry,  illustrated  in 
Figure  15,  consists  of  an  upper  and  lower  wall,  each  diverging 
at  15  degrees.  The  initial  plane  is  at  an  axial  distance  of 
one  meter,  and  the  final  plane  is  at  four  meters.  The  walls 
diverged  at  15  degrees  from  the  centerline.  For  an 
explanation  of  the  exact  solution  of  the  source  flow,  see 
Appendix  0  in  reference  3. 

The  normalized  static  pressure  error  was  computed  at  each 
nodal  location  as  follows: 

normalized  static  pressure  error  =  (75) 

^exact 


Interior  Point  Error  Convergence.  The  interior  point 
error  convergence  was  accomplished  by  analyzing  the  flow  along 
a  single  grid  point  location  that  corresponded  to  the 
centerline  of  the  source  flow  for  methodll,  method21,  and 
method22.  For  each  method,  a  configuration  of  11  nodes,  21 
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nodes,  and  31  nodes  at  each  computational  plane  was  studied. 
The  absolute  value  of  the  normalized  static  pressure  error  was 
integrated  using  trapezoidal  integration.  The  results  can  be 
seen  in  Table  II. 


l.ible  II.  Integrated  percent  error  in  static  pressure  along 
centerline . 


methodll 

method21 

method22 

11  nodes 

0.847 

0.401 

0.531 

21  nodes 

0.394 

0.171 

0.247 

31  nodes 

0.234 

0.117 

0.168 

For  methodll,  see  Figure  16,  increasing  the  number  of 
nodes  from  11  to  21  results  iu  a  53  percent  reduction  in 
integrated  error.  This  is  consistent  with  a  first-order, 
finite-difference  method. 

Figure  17  shows  the  error  convergence  for  method  21. 
Again,  increasing  the  number  of  nodes  from  11  to  21  results  in 
a  57  percent  reduction  in  error.  While  a  57  percent  error 
reduction  is  slightly  better  than  a  first-order,  finite- 
difference  convergence,  it  falls  short  of  the  75  percent 
reduction  that  would  be  achieved  for  a  second-order,  finite- 
differenced  equation. 

The  error  convergence  for  method22  can  be  seen  in 
Figure  18.  Increasing  the  number  of  nodes  from  11  to  21 
results  in  an  integrated  error  reduction  of  53  percent. 
Interestingly,  this  is  the  same  reduction  seen  in  methodll. 


and  is  very  close  to  raethod22's  value.  Again,  the  error 
reduction  falls  short  of  the  75  percent  value  that  would  be 
achieved  by  a  second-order  finite  differenced  equation. 

As  stated  before,  the  FDS  method  is  similar  to  finite- 
differencing  up  to  the  point  where  the  flux  differences  are 
split.  After  the  flux  differences  are  split,  only  a  portion 
of  the  information  from  a  particular  Riemann  location  might  be 
used  to  compute  new  values  at  the  next  computational  plane. 
A  finite-differencing  scheme  would  use  all  the  information, 
depending  on  the  stencil,  in  determining  the  values  at  the 
next  plane.  Because  of  this  splitting,  and  using  different 
amounts  of  information  from  the  same  location,  FDS  does  not 
directly  relate  to  finite-differencing. 

Boundary  Point  Error  Convergence.  The  boundary  point 
error  convergence  was  done  by  looking  at  the  flow  along  the 
upper  wall  of  the  source  flow  for  methodll,  method21,  and 
method22.  For  each  method,  a  configuration  of  11  nodes,  21 
nodes,  and  31  nodes  at  each  computational  plane  was  studied. 
The  results  can  be  seen  in  Table  III. 


Table  III.  Integrated  percent  error  in  static  pressure  along 
upper  boundary. _ 


methodll 

method21 

method22 

11 

nodes 

3.554 

1.658 

0.629 

21 

nodes 

1.587 

0.771 

0.410 

31 

nodes 

1.008 

0.482 

0.348 
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For  methodll,  see  Figure  19,  increasing  the  number  of 
nodes  from  11  to  21  results  in  a  55  percent  reduction  in 
error . 

Figure  20  shows  the  error  convergence  for  method21. 
Again,  increasing  the  number  of  nodes  from  11  to  21  results  in 
a  53  percent  reduction  in  error.  This  is  approximately  the 
same  error  convergence  as  methodll  above. 

The  error  convergence  for  method22  can  be  seen  in 
Figure  21.  Increasing  the  number  of  nodes  from  11  to  21 
results  in  an  integrated  error  reduction  of  only  35  percent. 
Method22  on  the  boundary  is  the  only  method  that  shows  the 
true  error  convergence  for  the  second-order  boundary  since  the 
boundary  effects  for  method22  at  the  centerline  are 
diminished.  This  lesser  reduction  of  error  at  the  boundary 
could  be  the  result  of  extrapolating  a  split  flux  difference 
to  a  location  outside  the  boundary,  and  weighing  it  evenly 
with  split  flux  difference  information  coming  from  a  Riemann 
location  just  inside  the  boundary.  In  effect,  nonphysical 
information  is  being  used  equally  with  physical  information 
near  the  boundary. 

In  summary  of  both  interior  and  boundary  integrated  error 
comparison,  all  three  methods  appear  to  have  a  near  first- 
order  finite-difference  convergence  behavior.  The  expected 
75  percent  error  reduction  for  a  second-order  method  was  not 
demonstrated  for  the  FDS  method.  This  is  caused  by  the 
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splitting  of  the  fluxes,  and  using  the  information  differently 
than  a  finite-differencing  method  would  use  it.  As  a  result 
of  the  grid  refinement  study,  second-order  accuracy  for  the 
FDS  method  cannot  be  verified. 

Method  Comparison 

To  further  compare  the  three  methods,  flow  traveling 
through  an  oblique  shock  wave  was  studied.  The  oblique  shock 
study  is  a  more  demanding  and  significant  test  than  the  source 
flow  study.  For  the  shock  study,  a  uniform  flow  of  Mach 
number  2.2  enters  a  channel,  region  1,  with  initially  parallel 
walls,  see  Figure  22.  As  the  flow  proceeds  downstream,  it 
encounters  a  ramp  at  the  bottom  wall  with  a  deflection  angle 
of  10  degrees.  This  ramp  creates  an  incident  oblique  shock 
wave  that  turns  the  flow  parallel  to  the  ramp,  forming  a 
second  area  of  uniform  flow,  which  is  represented  by  region  2 
in  the  figure.  Finally,  the  flow  is  turned  back  parallel  to 
the  wall  by  an  oblique  shock  reflection.  This  results  in 
uniform  flow  in  region  3.  All  three  methods  were  analyzed 
using  51  nodes  on  each  computational  plane.  Table  IV  gives 
the  exact  conditions  for  the  three  regions: 
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Table  IV.  Exact  values  in  oblique  shock  study. 


Property 

Region  1 

Region  2 

Region  3 

static 

pressure 

(Pa) 

206,842.0 

344,829.8 

548,689.3 

static 

temperature 

(K) 

1500.0 

1666.8 

1833.4 

Mach  number 

2.2 

1.885555 

1.582704 

specific 
heat  ratio 

1.25 

1.25 

1.25 

gas  constant 
(J/Kg/K) 

332.56 

332.56 

332.56 

Interior  Point  Method  Comparison.  A  graph  showing  plots 
of  all  three  methods  computed  across  an  oblique  shock  wave, 
along  the  streamwise  grid  location  of  j=41,  can  be  seen  in 
Figure  23.  At  first  glance,  it  appears  that  all  three  methods 
predicted  identical  pressures  for  a  given  value  of  x.  What 
this  means  is  that  if  the  pressures  were  calculated  at  every 
value  of  X,  the  values  computed  for  all  three  methods  would 
appear  on  the  same  curve.  But,  in  a  numerical  scheme,  values 
are  only  computed  at  discrete  locations,  which  is  determined 
by  the  step  size.  Both  second-order  interior  methods  compute 
the  exact  same  value  with  identical  step  sizes,  but  the  first- 
order  interior  method  computes  a  different  step  size. 

The  reason  that  all  three  methods  lie  very  near  the  same 
curve  is  because  the  magnitude  of  the  second-order  correction 
term  is  very  small  as  compared  to  the  first-order  solution 
that  it  is  added  to.  To  explain  why  the  second-order 


48 


correction  is  small  compared  to  the  first-order  solution,  a 
discussion  of  the  ideal  Riemann  problem  for  the  oblique  shock 
wave  is  in  order.  For  an  idealized  Riemann  problem,  on  each 
computational  plane  the  shock  wave  is  captured  by  two  nodes 
and  only  passes  through  one  Riemann  location  that  is  between 
the  nodes.  See  Figure  24.  In  this  case,  the  Riemann  location 
above  j,  at  j+1/2,  contains  nonzero  values  because  it 
separates  two  different  uniform  flow  regions,  at  j  and  j+1. 
Eg  (38)  can  be  used  to  demonstrate  this: 

(dEl)  ~  [  wave  3  wave  2  wave  1  ^  j+1/2 


As  can  be  seen  in  Eq  (38),  which  represents  the  total 
difference  between  nodes  j  and  j+1,  the  flux  differences,  at 
j+1/2,  would  contain  nonzero  values.  A  similar  argument  can 
be  made  about  the  Riemann  location  that  would  be  below  j,  at 
j-1/2.  Because  the  uniform  flowfields  at  j  and  j-1  are  equal, 
there  is  no  flux  difference  at  j-1/2.  Thus,  the  Riemann 
location  at  j-1/2,  is  in  uniform  flow  and  contains  only  zero 
values.  The  second-order  correction  computes  the  difference 
between  the  positive  biased  split  flux  differences  above  and 
below  the  node,  and  the  differences  between  the  negative 
biased  split  flux  differences  above  and  below.  The 
differences  of  the  positive  and  negative  split  flux 
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differences  from  above  and  below  would  give  values  equal  in 
magnitude  to  those  at  j+1/2  since  the  values  of  the  split  flux 
differences  below  are  equal  to  zero.  These  differences  of 
split  flux  differences  are  multiplied  by  a  coefficient  for  the 
second-order  corrections  and  are  then  added  to  the  first-order 
solutions.  In  the  idealized  case,  values  on  the  order  of  the 
coefficients  would  be  added.  In  the  case  where  the  shock  is 
smeared  over  four  or  five  nodes,  the  values  at  the  Riemann 
locations  are  also  smeared.  This  means  that  the  values  in  two 
adjacent  Riemann  locations  are  very  near  one  another.  As  a 
result,  the  differences  of  the  positive  and  negative  split 
flux  differences  above  and  below  the  node  are  very  small,  and 
when  multiplied  by  the  coefficients,  they  become  much  smaller 
in  comparison  with  the  first-order  solution.  Since  the 
second-order  corrections  that  are  added  to  the  first-order 
solution  is  much  smaller  than  the  first-order  solution,  all 
three  methods  appear  to  predict  values  along  the  same  curve. 

Boundary  Point  Method  Comparison.  A  graph  showing  plots 
of  all  three  methods  computed  on  the  upper  wall  across  a  shock 
reflection  can  be  seen  in  Figure  25.  Methodll  and  method21 
show  a  monotonic  convergence  behavior  across  the  shock.  This 
is  the  expected  result  since  both  methods  are  first-order 
accurate  on  the  boundary.  As  noted  above  in  the  interior 
point  method  comparison,  these  two  methods  appear  to  predict 
values  that  fall  almost  on  the  same  curve  with  the  only 
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difference  being  step  size.  For  iiiethod22,  the  solution  moves 
in  the  wrong  direction  just  before  the  shock,  and  slightly 
overshoots  the  exact  solution  after  the  shock.  The  initial 
drop  in  the  value  before  the  shock  is  caused  by  the  linear 
extrapolation  of  the  values  that  are  located  at  j-1/2  and 
j-3/2.  Recall  that  this  linear  extrapolation  was  needed  to 
get  values  for  split  flux  differences  at  a  fictitious  Riemann 
location  just  outside  the  boundary.  Information  was  needed 
there  to  compute  the  differences  of  split  flux  differences  for 
the  second-order  corrections.  The  problem  arises  because  the 
Riemann  location  at  j-3/2  may  be  in  region  2,  behind  the  shock 
wave,  but  j-1/2  may  still  be  influenced  by  region  1,  as 
illustrated  in  Figure  26.  This  results  in  erroneous 

information  being  extrapolated  to  the  fictitious  Riemann 
location.  The  same  problem  arises  when  the  computational 
plane  is  just  behind  the  shock  reflection.  Riemann  location 
j-1/2  is  influenced  by  region  3,  but  location  j-3/2  is  still 
in  region  2.  As  a  result,  the  second-order  boundary 

correction  looses  accuracy  around  a  shock  reflection. 


Figure  15.  Geometry  for  planar  supersonic  source  flow  (3:237) 
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Figure  16.  Source  flow  centerline  error  comparison  for  methodll 
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Figure  17.  Source  flow  centerline  error  comparison  for  method21 
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Figure  18.  Source  flow  centerline  error  comparison  for  method22. 
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Figure  19.  Source  flow  upper  boundary  error  comparison  for  raethodll 
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Source  flow  upper  boundary  error  comparison  for  method21 
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Figure  21.  Source  flow  upper  boundary  error  comparison  for  method22 


22.  Geometry  for  shock  wave  reflection  study  (3:44). 
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Figure  23.  Method  comparison  for  flow  across  an  oblique  shock  wave  at  an  interior 
point . 


X 

Figure  24.  Riemann  problem  for  a  shock  wave  (3:180) 
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Figure  26.  Linear  extrapolation  across  a  shock  at  a  boundary  (3:44). 
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VII .  Summary 


Second-order  correction  terms  were  added  to  an  existing 
computer  program  that  is  based  on  a  first-order  accurate  flux- 
difference-splitting  method.  Individual  analyses  were 
performed  in  an  effort  to  verify  the  second-order  accuracy  of 
the  solution  and  to  determine  the  benefits  associated  with  the 
second-order  method. 

A  grid  refinement  study  was  performed  in  an  attempt  to 
verify  a  second-order  accuracy  in  the  solution  obtained  by 
adding  the  second-order  correction  terms  to  the  first-order 
solution.  Second-order  accuracy  could  not  be  confirmed. 
Additionally,  it  was  noted  for  all  three  methods  evaluated 
that  doubling  the  number  of  nodes  resulted  in  an  error 
reduction  of  approximately  50  percent.  The  only  exception  to 
this  50  percent  reduction  was  seen  in  the  method22  analysis  at 
the  boundary,  which  had  a  reduction  of  only  35  percent.  This 
boundary  behavior  for  method22  resulted  from  extrapolating 
information  from  two  interior  Riemann  locations  that  were  in 
different  flow  regions  to  a  fictitious  Riemann  location 
outside  the  boundary,  and  then  weighing  the  fictitious 
information  equally  with  real  information  from  the  location 
just  inside  the  boundary. 

During  the  oblique  shock  reflection  study  it  was  observed 
that,  as  a  result  of  shock  smearing,  the  second-order 
correction  for  the  interior  points  was  very  small  in  magnitude 
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with  respect  to  the  first-order  solution.  Consequently,  all 
three  methods  predicted  values  that  appeared  to  lie  on  the 
same  curve.  The  only  deviation  was  that  the  step  size  for  the 
second-order  methods  is  different  than  that  of  the  first-order 
method.  The  shock  reflection  study  on  the  boundary 
demonstrated  that  methodll  and  method21  predicted  similar 
values.  This  was  expected  since  the  two  methods  are  first- 
order  accurate  on  the  boundary.  Again,  the  only  deviation  was 
in  the  step  size  and  the  axial  location  of  the  computational 
planes.  Method22  performed  poorly  on  the  boundary.  The 
reason  for  the  poor  performance  for  method22  was  because  it 
uses  information  extrapolated  from  the  first  two  interior 
Riemann  locations  to  assign  information  to  a  fictitious 
Riemann  location  outside  the  boundary.  The  problem  lies  in 
the  fact  that  the  method  extrapolates  across  the  oblique  shock 
wave,  resulting  in  erroneous  values  at  the  nonphysical  Riemann 
location.  This  anomaly  occurs  both  in  front  of  the  shock 
reflection  and  behind  it. 

As  a  result  of  the  above  observations  and  conclusions, 
and  with  the  knowledge  that  the  first-order  method  has  a 
faster  computation  time  and  is  as  accurate  as  many  second- 
order  finite  differencing  schemes,  Taylor  (8:108),  the 
recommendation  is  to  perform  future  analysis  using  the  first- 
order  method. 
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Appendix  A;  Weak  Conservation  Law  Form 
of  the  Governing  Equations 


The  weak  conservation  form  uses  metrics  at  the  Riemann 
locations  and  is,  therefore,  consistent  with  the  FDS.  This 
translation  of  the  metrics  is  a  logical  step  because  the  split 
flux  differences,  which  pass  information  to  plane  i+1,  are 
computed  at  the  Riemann  locations.  Converting  the  governing 
equations  into  weak  conservation  form  is  done  by  moving  the 
metrics  in  the  governing  equations,  Eq  (10),  into  the 
derivative,  and  differencing  the  flux/metric  product, 
Hindman  (6:113).  A  mathematical  manipulation  is  require  to 
move  the  metrics  into  the  differentiation  that  occurs  in 
Eq  (10),  which  is  repeated  here  for  convenience: 


d{S)  _  „  d{S)  din 

ac  ^  ail  ari 


(10) 


Using  the  product  rule  of  partial  differentiation,  it  can  be 
demonstrated  that  the  metrics  can  be  moved  direv-wly  into  the 
differentials  on  the  right-hand  side  of  Eq  (10)  as  follows: 


a  , 


«)  =  n; 


ari 


(76) 


Manipulating  the  last  term  in  Eq  (76)  shows 
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(77) 


*  =  ^XT,  «  =  =  (D^S  =  0 


Substituting  Eg  (77)  into  Eg  (76)  and  reversing  the  terms 


(78) 


Eguation  (78)  corresponds  to  the  first  term  on  the  right-hand 
side  of  Eg  (10).  A  similar  manipulation  can  be  performed  on 
the  last  term  in  Eg  (10).  Therefore,  it  has  been  demonstrated 
that  the  metrics  can  be  moved  into  the  differentials  while 
retaining  the  integrity  of  the  governing  eguations.  Thus, 
through  a  similar  derivation  of  Eg  (59),  Eg  (10)  can  be 
transformed  into  the  first-order  flux-difference-splitting, 
weak  conservation  law  form: 


nr  =  [cix  -  (’i. 

-  A(  [(,,  dr-),.,,,  -  (I,  dr-) 


(79) 
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